% Taking the firm-side calibration as given, this script does the following:
% 1) For a given markup at Rz = 5% targeted by the policymaker, it finds a value of
% Hotelling parameter xi that would rationalize it. 
% 2) It finds a new equilibrium by checking which contract on the modified
% iso-profit maximises agent's perceived utility. 
% 3) It calculates the impact on firm's profits and consumer welfare,
% relative to laissez-faire. 



clear
load('recalibrated_data_with_eval.mat')

fileID = fopen('Policy_2_competition_2021.txt','wt')


% Specifiy markup after the intervention: 
markup = 0.30;

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Policy experiment with a NEW MARKUP equal to %f  at Rz = 0.05 \r\n'],[markup]);


rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;
age_ = mydata.age;

THETAtarget = 0.05;
totalcosttarget = 0.005;
adminratiotarget = 0.50;
%markuptarget = 0.2;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_naif_target = mydata.(suffix);
suffix = strcat('fmax_',num2str(THETAtarget*10000));
fmax_target = mydata.(suffix);


% Set the calibrated values:
gamma1 = 0.5;
c1 = adminratiotarget*totalcosttarget*mean(avgZ_naif_target) / mean(avgZ_naif_target.^(gamma1));
gamma2 = 5.75;
c2 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2);
%phi = 0.35;

% And for later comparisons with laissez-faire:
lf_profit = 15211.219049;
lf_actual_utility = 654.673018;


% First, calculate phi associated with the new markup
costs = c1*avgZ_naif_target.^(gamma1) + ones(1,71)*c2*(1+THETAtarget-R)^(gamma2);
phi = 0;
fee = phi*fmax_target;
while (fee/mean(costs) - 1) < markup
    phi = phi + 0.005;
    fee = phi*fmax_target;
end

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['The requested markup of %f is associated with phi of %f \r\n'],[markup,phi]);


% Second, back out new xi (which the intervention effectively targets)
xifound = 0;

% 1a) Calculate utility from the (supposed) equilibrium contract
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETAtarget*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAtarget*10000),'_TC');
avgTC_TC = mydata.(suffix);
    
bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

avgTC_fee = avgTC_TC - phi*fmax_target*ones(1,71);
    
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
    
lifetime_equilibrium = instant(1);
for t = 2:71
    lifetime_equilibrium = lifetime_equilibrium + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

% 1b) Check that the firms don't want to deviate to lower/higher fees for target Rz

gridlength = 11;
gridjump = 0.05;

ProfitPVvec = zeros(1,gridlength);
expected_profits = zeros(1,gridlength);
xi = 0.0001;

for j = 1:gridlength
    
    % First, profit from an interaction
    fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*phi*fmax_target;
    
    THETA = THETAtarget;
    avgZ_ = avgZ_naif_target;
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

    profits = ones(1,71)*fee - costs;

    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
    
    ProfitPVvec(j) = ProfitPV;
    
    % Second, probability of interaction
    avgTC_fee = avgTC_TC - fee*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_fee = instant(1);
    for t = 2:71
        lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    probabilityraw = 0.5 + (lifetime_fee - lifetime_equilibrium)/(2*xi);
    probability = max(0,min(1,probabilityraw));
    
    expected_profits(j) = ProfitPV*probability;
end

if max(expected_profits) == expected_profits(ceil(gridlength/2))
    xifound = 1;
    fprintf(fileID, ['************** \r\n']);
    fprintf(fileID, ['For xi = %f, the calibrated fee of %f is Hotelling equilibrium \r\n'],[xi,phi*fmax_target]);
else
    xi = 1.25*xi;
end

counter = 1;
while xifound < 1 && counter < 1000
    expected_profits = zeros(1,gridlength);
    for j = 1:gridlength
        
        fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*phi*fmax_target;
    
        avgTC_fee = avgTC_TC - fee*ones(1,71);
    
        instant = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
            else
                instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        lifetime_fee = instant(1);
        for t = 2:71
            lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
        end
    
        probabilityraw = 0.5 + (lifetime_fee - lifetime_equilibrium)/(2*xi);
        probability = max(0,min(1,probabilityraw));
    
        expected_profits(j) = ProfitPVvec(j)*probability;
    end

    if max(expected_profits) == expected_profits(ceil(gridlength/2))
        xifound = 1;
        fprintf(fileID, ['**************\r\n']);
        fprintf(fileID, ['For xi = %f, the calibrated fee of %f is Hotelling equilibrium \r\n'],[xi,phi*fmax_target]);
    else
        xi = 1.25*xi;
    end
    
    counter = counter + 1;
end


% 1c) Profit level at the target
THETA = THETAtarget;
avgZ_ = avgZ_naif_target;
    
costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

profits = ones(1,71)*phi*fmax_target - costs;

maxprofit = profits(1);
for t = 2:71
    maxprofit = maxprofit + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
end


% 1d) Need to re-calculate the iso-fees
for THETA = mydata.THETA_grid
    if abs(THETA - THETAtarget) < eps
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = phi*fmax_target;
    else
        suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
        avgZ_ = mydata.(suffix);
                        
        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
                        
        fee = 100;
        profits = ones(1,71)*fee - costs;
        ProfitPV = profits(1);
        for t = 2:71
            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
        end
                    
        while maxprofit > ProfitPV
            fee = fee + 6;
            profits = ones(1,71)*fee - costs;
            ProfitPV = profits(1);
            for t = 2:71
                ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
            end
        end
                        
        suffix = strcat('fmax_',num2str(THETA*10000));
        fmax = mydata.(suffix);
                        
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = min(fee-3,fmax);
    end
end

                    
% 1e) Calculate profits for each Rz (just to double-check the iso-profit condition)
fprintf(fileID, ['**************\r\n']);
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    
    % This is just needed for calculating out-of-eqm outcomes
    %if abs(THETA - THETAtarget) < eps
        %costs500 = costs;
    %end
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
         
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
             
    fprintf(fileID, ['For Rz = %f, the PV of profits (on the iso-profit curve) is %f \r\n'],[THETA,ProfitPV]);
end


% 2) Given the new markup at Rz = 5% and recalculated iso-fees, which
% contract maximises perceived utility and is therefore offered in
% equilibrium?
fprintf(fileID, ['**************\r\n']);

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end
                
perceived_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
    avgZ_TC = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_TC');
    avgX_TC = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
    avgTC_TC = mydata.(suffix);
                    
    suffix = strcat('isofee_',num2str(THETA*10000));
    isofee = workspace.(suffix);
    avgTC_TC = avgTC_TC - ones(1,71)*isofee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_TC(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_TC(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    perceived_lifetime_iso_grid = [perceived_lifetime_iso_grid lifetime];
    
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the PERCEIVED lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
    if abs(THETA - THETAtarget) < eps
        lifetime_at_target = lifetime;
    end
end

[market_benchmark, id_market] = max(perceived_lifetime_iso_grid);
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Now, Rz =  %f maximises perceived utility and is offered in market equilibrium \r\n'],[mydata.THETA_grid(id_market)]);



% 2) If Rz offered in new equilibrium differs from laissez-faire, need to
% double-check no profitable price deviations condition:
% There should be no incentive to price-deviate from
% iso_fee(id_max_profit) by +/- 25%: 
THETA = mydata.THETA_grid(id_market);

suffix = strcat('isofee_',num2str(THETA*10000));
fee_supp_eq = workspace.(suffix);
suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
avgZ_ = mydata.(suffix);

suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETA*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
avgTC_TC = mydata.(suffix);

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

% Calculate perceived utility at current eqm candidate
avgTC_fee = avgTC_TC - fee_supp_eq*ones(1,71);    
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
lifetime_new_eq = instant(1);
for t = 2:71
    lifetime_new_eq = lifetime_new_eq + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

% Verify no profitable price deviation: 
gridlength = 11;
gridjump = 0.05;

ProfitPVvec = zeros(1,gridlength);
expected_profits = zeros(1,gridlength);

for j = 1:gridlength
    % First, profit from an interaction
    fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*fee_supp_eq;
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
    ProfitPVvec(j) = ProfitPV;
    
    % Second, probability of interaction
    avgTC_fee = avgTC_TC - fee*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_fee = instant(1);
    for t = 2:71
        lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    probabilityraw = 0.5 + (lifetime_fee - lifetime_new_eq)/(2*xi);
    probability = max(0,min(1,probabilityraw));
    
    expected_profits(j) = ProfitPV*probability;
end

[max_profit,adjustment_id] = max(expected_profits);
optimal_adjustment = (1 - floor(gridlength/2)*gridjump + (adjustment_id-1)*gridjump);
fee_eq = optimal_adjustment*fee_supp_eq;
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['For xi = %f and Rz = %f, there is no profitable price-based deviation for fee = %f  \r\n '],[xi,THETA,fee_eq]);
fprintf(fileID, ['(Adjustment (if any) with optimal multiplier of  %f implemented)  \r\n '],[optimal_adjustment]);

% Note: Need to multiply max_profit by 1/0.5 to correct for eqm
% probability of interaction:
max_profit = max_profit*2;




% 3) Impact on firms profits
profitchange = (maxprofit - lf_profit)/lf_profit;

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Relative to laissez faire, the change in firms profits is %f percent\r\n'],[100*profitchange]);



% 4) Calculate the ACTUAL consumer welfare assocviated with the new
% equilibrium contract 

fprintf(fileID, ['**************\r\n']);
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

actual_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid(id_market)
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_naif');
    avgX_ = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
    avgTC_ = mydata.(suffix);
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
                    
    avgTC_ = avgTC_ - ones(1,71)*fee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    %actual_lifetime_iso_grid = [actual_lifetime_iso_grid lifetime];
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the ACTUAL lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
end


% What is the CE welfare loss?
utility_multiplied = lifetime;
efficient = lf_actual_utility;
multiplier = 1.00;

if utility_multiplied < efficient
    while utility_multiplied < efficient 
        instant_multiplied = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
            else
                instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        utility_multiplied = instant_multiplied(1);
        for t = 2:71
            utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
        end
        multiplier = multiplier + 0.0001;
    end

    fprintf(fileID, ['The CE welfare LOSS relative to laissez faire is %f percent\r\n'],[100*(multiplier-1.0001)]);

else
    while utility_multiplied > efficient
        instant_multiplied = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
            else
                instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
        
        utility_multiplied = instant_multiplied(1);
        for t = 2:71
            utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
        end
        
        multiplier = multiplier - 0.0001;
        
    end
    fprintf(fileID, ['The CE welfare GAIN relative to laissez faire is %f percent\r\n'],[100*(1-0.0001-multiplier)]);
end

        
    



    
    
                    
                    
